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Abstract. We propose a new dynamical system formalism for the analysis of f(R) cosmolo¬ 
gies. The new approach eliminates the need for cumbersome inversions to close the dynamical 
system and allows the analysis of the phase space of f(R )-gravity models which cannot be 
investigated using the standard technique. Differently form previously proposed similar tech¬ 
niques, the new method is constructed in such a way to associate to the fixed points scale 
factors, which contain four integration constants (i.e. solutions of fourth order differential 
equations). In this way a new light is shed on the physical meaning of the fixed points. 
We apply this technique to some f(R) Lagrangians relevant for inflationary and dark energy 
models. 
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1 Introduction 

Since the first formulation of General Relativity (GR), many extensions of the original Ein¬ 
stein equations have been investigated. The reasons of the interest in such theories are quite 
disparate: from the first attempts to unify geometrically the electromagnetic and gravitational 
interaction started by Weyl [1], to the understanding of the corrections to the gravitational ac¬ 
tion typical of quantum held theory in curved spacetime and fundamental unification schemes 
[2], to the attempts to give a complete geometric explanation of the dark phenomenology. 

Among these extensions, the class of theories called f(R )-gravity [3] is the simplest 
realisation of higher order gravity (order four) and has an important role as a natural model 
for inflation [4], More recently, it was shown to have also an interesting (and very debated) 
role as geometrical dark energy model [8]. In dealing with these theories, the necessity to solve 
fourth order differential equations is the source of the difficulties in the true understanding of 
the features of their cosmology. Such problems were the origin of the development of a series 
of methods to indirectly analyse the physical properties of f(R) cosmologies. 

Dynamical System Approach (DSA) has proven to be one of the most effective of such 
methods. DSA has a number of different realisations [5, 6]. In the following we shall consider 
the one proposed in the 1970’s by Collins [7] and successively developed by Wainwright, Ellis 
and Uggla. This version of DSA is defined in terms of dimensionless, expansion normalised 
variables and allows a very clear physical interpretation of the results. DSA has had a key role 
in analysing Bianchi models [9] and minimally coupled scalar tensor [10] in the context of GR 
and has allowed the exploration of the cosmology of a number of modifications of Einstein’s 
theory [11, 12]. In [11], the model f{R) = xR n was analysed with this method for the fist 
time. Later on the method was extended in [13, 14] to the case of a generic f(R). DSA 


- 1 - 




allowed for the first time to analyse in detail the phase space for /(/?) cosmologies making 
some generals statements on these cosmologies and revealing a number of interesting attractor 
solutions. 

In spite of this success the method above has some unsatisfactory aspects. First of all 
the possibility to analyse a given f(R) theory depends on the exact resolution of an algebraic 
equation of generic order, or, in the most complicated cases, a transcendental equation. Such 
operation not only limits the set of possible Lagrangians which can be analysed with DSA, 
but also introduces a number of singularities, so that in many cases the dynamical system is 
not of class C 1 . 

Another important problem is that the choice of the dynamical system variables was 
made in such a way to obtain in the fixed points solutions with only two integration constants. 
This implies that these solutions correspond to the general solutions of the fourth order 
cosmological model where two integration constants have been set to zero. This is problematic 
because without knowledge of the full solution in a fixed point it is impossible to characterise 
the correct behaviour of these cosmologies when these fixed points are nodes. 

In addition, since the dynamical system variables are not always independent from each 
other, the dynamical system might present fixed points which can correspond to inconsistent 
conditions in the cosmological equations i.e. to have f(R) = 0 and R = 0 for a function / 
for which /(0) 0. The presence of these points must be a spurious effect due to the way in 
which he DSA is constructed. 

In this paper we propose a new DSA to deal with the cosmology of f(R) gravity. This 
method is built in such a way to avoid the necessity of solving exactly algebraic/transcendental 
equations and to give, in the fixed points, solutions of the cosmological equations which 
contain four integration constants. We will show that the new formulation contains the 
results of the old method but reveals unsuspected additional features of the evolution of f(R) 
cosmologies. 

The paper is organised as follow. In section 2 we will give the basic equations. In 
section 3 we will introduce the original DSA for f(R)- gravity. Section 4 is dedicated to the 
construction of the new DSA, and in Section 5 the new method is applied to some interesting 
model of f(R) gravity. Section 6 is dedicated to the conclusions. 

Unless otherwise specified, natural units {K = c = ks = 87 tG = 1) will be used through¬ 
out this paper, Latin indices run from 0 to 3. The symbol V represents the usual covariant 
derivative and d corresponds to partial differentiation. We use the (—, +, +, +) signature and 
the Riemann tensor is defined by 

R\cd = w a bd>c - W\ c4 + W e bd W a ce - W f bc W a df , (1.1) 

where the W a bd are the Christoffel symbols (i.e. symmetric in the lower indices), defined by 

hb bd = ~^9 ( 9be,d T 9ed,b 9bd,e) • (1-2) 

The Ricci tensor is obtained by contracting the first and the third indices 

Rab = g Cd Racbd ■ (1-3) 

Finally the Hilbert-Einstein action in the presence of matter is given by 

A = J dx 4 ^g ^R + L m . (1.4) 
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2 Basic Equations 


In this paper we are going to deal only with metric f(R) theories, which are characterized by 
the Action 

S = j (fixyf^g [f(R,a,j3...) + C m \ , ( 2 . 1 ) 

where C m represents the matter contribution. In general the function / is considered an 
analytic function of the Ricci scalar R and contains a set of additional dimensional parameters 
indicated by barred Greek letters. 

Varying the action with respect to the metric gives the generalisation of the Einstein 
equations: 

f'G ab = T2 + l -g ab (/ - Rf) + ( g a c g b d - g ab g cd )V C V d f , (2.2) 

where G ab is the Einstein tensor, / = f(R , a, /?...), f = —— ’ ~^ —-, and TM = — ^ m ^ 

dR yj-g og ab 

represents the stress energy tensor of standard matter. These equations reduce to the standard 
Einstein held equations when f(R, a, (3...) = aR with a = 1/2. 

Our treatment will consider only homogeneous and isotropic spacetimes i.e. Friedmann 
Lemaitre Robertson Walker (FLRW) metrics: 


ds 2 


— dt 2 + a 2 (t.) 


dr 2 

1 — kr 2 


+ r 2 (d6 2 + sin 2 Odcj) 2 ) , 


(2.3) 


where a is the scale factor and k the spatial curvature. We also assume that the cosmic fluid 
is a prefect fluid with equation of state p = wg with 0 < w < 1. It is common to write the 
held equations (2.2) in the metric (2.3) as two equation resembling the Raychaudhuri and the 
Friedmann equations in GR 


« 2 + £ = ^7 { \ [f'R - /] - 3Hf' + , 

2rr + ff 2 + W-i|l [t'R -/]+/'- 3 Hf + p m \ , 


where 


R = 6 | H + 2H 2 + 4 


(2.5) 


H = d/a, the prime represents the derivative with respect to R and the “dot" is the derivative 
with respect to t. The two equations (2.4) are not independent: the second can be obtained 
deriving the hrst with respect to the cosmic time t once the Bianchi identities for T™ is 
considered. In FLRW these identities take the form 


Am T -iH(g m T Pm) — 0 j 

which is the same as the GR energy conservation equation for the cosmic huid. 
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3 The Original Dynamical Systems Approach for f(R) gravity (in brief). 


Using the equations above we can formulate the cosmic evolution in terms of dynamical sys¬ 
tems [13, 14] (beware of the differences in signature!). Introducing the general dimensionless 
variables : 


= ±_ R f 

X Hf ’ V 6H 2 ’ 6H 2 f' ’ 

ft = V m K = k 

3H 2 f” a 2 H 2 ’ 

where /i m represents the energy density of a perfect fluid that is present in the model. As 
customary, we also define the logarithmic (dimensionless) “time variable” N = In a. Note 
that in choosing this time variable we are assuming that we represent the phase space for 
H > 0 i.e. we are considering only expanding cosmologies 1 . In this variables the cosmological 
equations (2.4) are equivalent to the autonomous system: 


dx 

dN 

dy 

dN 

dz 

dN 

dD 

dN 

0 


x(D-z- 2) - 2x 2 + 2 (y - 2 z) + (1 - 3w)D, 

y[(T -2)x + 2n-2z + 2], 

z (2 — 3x — 2z + 2Sl)z + x yY, 

n (2D, — 3x — 2z — l-3w;), 

1 + K + x + z — y — D . 


(3.2) 


The quantity Y is defined as 


Y = 


Rf 


(3.3) 


Since Y is a function of R only, the problem of obtaining Y = T(x,y, z,D) is reduced to the 
problem of writing R = R(x, y , z, D). This can be achieved from the definitions (3.1): 


v = Rf 

Z f 


(3.4) 


Solving the above equation for R allows one to write R in terms of y and z and close the 
system (3.2). It is clear that the properties of (3.4) determine the possibility of closing (and 
therefore analysing) the system (3.2) as well as some of the properties of this system i.e. the 
differential structure. One example is f(R) = RP exp (qR) for which 


T = 



(3.5) 


In this case it is evident that the system is not (7(1) as the curve y 2 — pz 2 = 0 is singular. 
This fact can have serious repercussions on the properties of the flow [14]. 

1 The contracting case can also be considered using the time variable M = —N. Since by definition almost 

all the variables (3.1) are invariant under a change of sign of H the phase space for H < 0 will have a strict 
resemblance with the one with H > 0. The H < 0 part of the phase space can be important to analyse, for 
example bouncing scenarios. In the following we will not consider this part of the phase space focusing only 
on expanding cosmologies, leaving this analysis for a future work. 
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4 The new Dynamical Systems Approach. 


We will now start to construct the new approach. Before we define the dynamical system 
variables we will need, however, two preliminary steps. The first one will concern the form 
of the action and the second one will be the introduction of new (cosmic) parameters. These 
(re-)definitions will be the cornerstones of the new method. 

4.1 The form of the action. 

In dealing with dynamical system it is crucial to gain an understanding (and control) over 
the dimensional structure of the theory we are considering. The reason is that the number of 
dynamical system variables needed for DSA will also depend on the number of dimensional 
constants present in the theory. To construct the new method, therefore, we will rewrite 
the action (2.1) in a special form. In particular we will introduce a constant Rq such that 
the product RRo is dimensionless. In addition, we will also introduce some dimensionless 
parameter in the form of Greek letters which will represent the ratio between the coupling 
constant of the additional invariants in the theory and (a power of) Rq. In this way the (2.1) 
can be written as 

S = j d A x^g[f(RR 0 ,a...)+C m ] , (4.1) 

where / has the same properties of the one in (2.1) and Ro will be assumed non-negative. 
The main reason behind the formulation above is that in this way any f(R) action contains 
only one dimensional constant. Therefore, instead of defining a dynamical variable for each 
dimensional constant (a, (3 ,...), we only need one dynamical variable related to Rq to analyse 
the phase space of actions of any complexity. In addition, this setting prevents the appearance 
of fixed points not consistent with the cosmological equations which are typical of the original 
DSA. We will use this formulation (4.1) of the action as a starting point in the construction 
of our new DSA. 

4.2 New cosmic parameters. 

Looking for a different way to constrain the properties of the cosmic fluids, Visser proposed 
a set of cosmic (or cosmographic) parameters [15-18] 

q = --H~ 2 , 8=—H-\ (4.2) 

a a a 

with which one is able to characterize completely a cosmological model. These quantities are 
directly related with the Taylor development of the scale factor and this property determines 
also our capability to measure them. With few exceptions [16], in the case of GR only the 
lowest order cosmic parameter have been so far fully exploited, but in the case of /(-R)-gravity 
the situation is different: higher order parameters become crucial and can be used to char¬ 
acterize the evolution of many important cosmological phenomena e.g. structure formation 

[19]- 

However, after a quick look to the f(R) cosmological equations written in terms of these 
quantities (e.g. [19, 20]), one soon realizes that this type of cosmographic parameters are 
not always the ideal objects to work with. The same happens when one tries to use them to 
formulate a dynamical system approach: although in principle the (4.2) are the ideal objects 
to construct the DSA they do not constitute always an advantageous set of variables. 
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For this reason, it is necessary to look for new sets of parameters which share the 
structure of the original cosmographic parameters (4.2) but, at the same time, are more 
suitable to deal with our specific problem. One possibility is to use the Hubble parameter to 
define the variables: 


H __ HH _ _ HH 2 
q ~ ir 2 '- 


(4.3) 


The variables above have been used in [21] to propose new ways to perform the reconstruction 
of exact cosmological solutions. Unfortunately also the (4.3) do not prove useful to implement 
a more powerful dynamical system approach. Another possibility, which will be adopted in 
the following, is to define 


H H H 2 

q ~lP ’ 




HH 

~W 


(4.4) 


This choice, apparently cumbersome, appears much simpler in terms of the logarithmic time 
N: 

H n . Hnn H nnn 

q = lT’ , = Hr' * = pT' <45) 

The (4.5) will be the second cornerstone of the mode we intend to propose. 

In terms of q,j,s the Ricci scalar and its derivatives read 


R = 6 


(q + 2) H~ + 


2 


R = 6tfj [j + q(q + 4)]i^--2 j, 

R = 6H 2 | [s + 4j(q + 1) + (q + 8)q“] H 2 + 2(2 — c l )^2 | > 


(4.6) 

(4.7) 

(4.8) 


and the cosmological equations (2.4) can be written as 

r2 f n j q H 4jn 


^ (1 + q) + 3 r + 12 ar*^- 


e r r 


2 iff 


+ 


+ 


r 

3 6i7 2 / (3) 

7 


[j + q(q+4)] =0, 


k 


377'[s + j(4q + 5) + q (q + 9) + 4q] + 6(1 — q)— / + 


+ q 2 + 4q) + 9 . + 2+ 2 _ 20 

a z 2 a* 


(4.9) 


respectively. Note that now the cosmological equations are equations for j and s instead of 
H and H. In fact, in (4.9), this last quantity has been substituted according to (4.4). 

In the system (4.9) the Ricci scalar only remains present in the function /. This means 
that the equations above should be supplied with an additional constraint given by the (4.6). 
This relation will allow to simplify considerably the final dynamical system. 


4.3 The General Method 

The first step to obtain and autonomous system of first order differential equations corre¬ 
sponding to the (4.9) is the definition of the dynamical variables. In the case of single fluid 


6 

















cosmology 2 , we choose the set of variables: 


R = 


R 

6W 2 ’ 


K = 


k 


J = 


l 

4’ 



Q = _ H _ 

' 3 H 2 f' 

A = R 0 H 2 . 


(4.10) 


Note that the variable associated to matter does not coincide exactly with the matter density 
parameter. This is a manifestation of the non-minimal coupling of matter and gravitation 
typical of /(ii)-gravity. Also, because of our definition of Ro, the variable A will be always 
non-negative. In addition to the above variables, we will introduce, like in the original DSA, 
the logarithmic time variable N. 

With this choice the cosmological equations (4.9) are completely equivalent to the au¬ 
tonomous system: 


^ = ^Q(Q - 3M + 9) - 2R + 4JJ + 4, 

^ = ^{[9M - 2(Q 2 + 9Q + 9J + 9)]Y - 6[4Q + 9(u> + 1)]}, 

^ + 9) + 9(1 + 1)]} 2 + 

+ + 4Q — 6R + 311(1 + 3 w) + 6] + 

O 1 

+ -^{9(3 - 2Q)R + 2(3 + 2Q)[Q(Q + 15) + 9(5J - 1)]}, 

54 

dQ 2Q 2 


dN 

dlC 

dN 

dA 

dN 



(4.11) 


together with the two constraints 


1 = n-K + R-X- 
2 

R = K + -Q + 2. 

3 



+ jr +1 


(4.12) 

(4.13) 


The first corresponds to the Hamiltonian (Friedmann) constraint, which guarantees the con¬ 
servation of matter energy, and the second is simply the definition of the Ricci scalar. The 
functions X = X (A, R), Y = Y (A, R), Z = Z (A, R) are defined respectively as 


X (A, R) = 


/ (R, A, a,...) 


Y(A,R) = 


24H 2 f" (R, A, a,...) 


6H 2 f' (R, A, a ,...) ’ v ’ ' /' (R, A, a,...) 5 

4 fill r™ a \ l 4 ' 14 ! 


Z (A, R) — 


96tf 4 /'"(R,A,a,...) 
/' (R, A, a,...) 


2 The approach can be trivially generalised to the multi-fluid case by adding a suitable number of fl variables 

each corresponding to the energy density of the given fluids. Such generalization does not add anything to 

the understanding of the method and it will not be pursued here. 
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These quantities represent the part of the system which depends on the form of the Lagrangian 
(4.1). Note that, differently from the approaches presented in [13, 14], in this version of the 
Dynamical Systems Approach no resolution of algebraic equation is required to close the 
system. This means that with the technique proposed here we can in principle analyze all 
f(R) models. 

The two constraints allow us to eliminate two variables (J and Q) and reduce the total 
system to: 


^ = 2R(K - R + 2) - ^(X + K - R - + 1), 

^ = D(2 - 3w + X + 3K - 3R - fi), 

H = 2K(K-R + 1), 
d a 

— = -2A(2 + K-R). 


(4.15) 


It is clear that this system posses a minimum of three invariant submanifolds: (i) = 0, (ii) 

K = 0 and (iii) A = 0. Depending on the form of Y one can also have a fourth invariant 
submanifold in M = 0. These invariant submanifolds can be imagined as “parts" of the phase 
space which have the property that any orbit that starts in them is trapped. They can have 
a very specific physical meaning. For example, the existence of the invariant submanifold 
K = 0 tells us that if we start in an orbit with zero spatial curvature we cannot evolve 
towards a positive or negative curvature parts of the phase space and that the existence of 
D = 0 implies that a vacuum cosmology remains vacuum (i.e. standard matter cannot be 
created or destroyed). The submanifold R = 0 represents, instead, the case in which the Ricci 
scalar is identically zero. Since the equation for R contains terms with Y at denominator, 
this submanifold can be singular. Finally, the submanifold A = 0 corresponds to the case in 
which the constant Rq is zero and it is of more difficult physical interpretation. It can be 
thought to correspond to the case in which the gravitational part of the action is identically 
zero i.e. there is no gravitational interaction 3 . Since A appears also in Y, the submanifold 
A = 0 can also be singular. 

The singular character of the invariant submanifolds R = 0 and A = 0, however, does 
not imply the absence of fixed points. The existence of fixed points that belong to a singular 
manifold is one of the exotic properties of dynamical systems which are not of order C(l). 
The coordinates of (and therefore the solution associated with) these points can be explained 
in terms of the limiting form of / for R, Rq —> 0. In fact, close to the fixed points in the 
R = 0 and A = 0 manifolds the phase space of a given theory / and the one of its limit for 
R,Ro —> 0 can be considered isomorphic and will admit the same fixed points. For these 
fixed points the stability analysis will not necessarily be the standard one. For example, some 
of the eigenvalues of these points might diverge. We will impose the condition that actual 
fixed points will exist if the equations at least of class C( 2) i.e. that the dynamical system 
equations, the cosmological equation and the eigenvalues will be finite for a given fixed point. 

3 At this point one might think that reformulating the dynamical system approach distinguishing the 
dimensional constant in front of the Hilbert-Einstein term and the one(s) of the higher order invariant(s) 
might ease the interpretation of this kind of invariant submanifolds. Indeed this can be done, but it does not 
add anything to the phase space analysis. For this reason we will rather keep using the approach presented 
above, which is more compact. 


- 8 - 



With the above setting, we are ready to analyse the phase space for a given form of 
/(R) 4 5 (and therefore a form of the functions X, Y, Z). The analysis can be done using 
the classical dynamical system theory i.e. obtaining the fixed points imposing that the N- 
derivatives of the dynamical variables are zero and using the Hartman-Grobman theorem to 
determine their stability’ 3 . Naturally, since the phase space is not compact one will need to 
perform an asymptotic analysis in order to obtain a full description of the dynamics. Such 
task, however, will be left for a future work. 

The solutions associated to the fixed points can be derived writing the modified Ray- 
chaudhuri equation in a fixed point, 


s = 


1 d3R _ 4 
HdNZ ~ 27 


4 

- 3 J* 


2Z* 


+ 33) Q* — 9M* + 27] + 
+ 9)Q* + 9(2 — M*)] + § 


+ 15K 


2z g 

27 Y ^(^* + ^(2 — -^*)] 2 + y ~ 


(4.16) 


+ D] + 


8 X* 

y7 


24J*Z 


40 * 


vy - + *^—(1 + 3u>) + 2(2 — M*), 

X * X * 


where the asterisk indicates the value of a variable in a fixed point 6 . In spite of being a 
third order differential equation the (4.16) is always solvable exactly for H. In addition, since 
(4.16) is a fourth order equation in a, the solution associated to the fixed points solution 
contains four integration constants (and up to three different “regimes”). As we will see, in 
the case in which a theory can be treated with both the new and the original method, there 
is correspondence between the fixed points. This suggests that the original method implicitly 
sets to zero some integration constants. We will discover that this can hide information on 
the actual meaning of the fixed point, particularly when it is a sink or a source. The other 
key cosmological quantities can be deduced by the cosmological equations once the (4.16) is 
solved. Specifically one has, from (2.6), 


^ = a~^ l+w \ 


(4.17) 


In literature one often defines the barotropic factor of the high order corrections considered 
as an effective fluid. 

Unfortunately, it will not always be possible to integrate the (4.16) exactly up to an 
expression for the scale factor. In this case we will rely on considerations on the structure of 
the equation for a and numerical integrations to obtain the behaviour of the scale factor and 
the other key quantities. For this reason we will limit ourselves to give the solutions for the 
scale factor in terms of plots. 

4 We should stress here that some general conclusions on the system (4.15) could be drawn without speci¬ 
fying the form of /. However this can be dangerous. Consider for example a case in which the function Y is 
proportional to R _1 . In this case the system would admit a fixed point ]R = 0,fi = 0,K = 0, A = 0 which 
would be not necessarily present in other cases. For this reason we will refrain form such speculations and 
will work only on a given forma on /. 

5 It is clear that, as some free parameters enter in the dynamical equations there will be values of these 
parameters for which the eigenvalues of the fixed points are zero. Such phenomena are called bifurcations (see 
e.g. [29]) and we will not explore them here. 

°It is important to stress here that the (4.16) represents the form that the Raychaudhuri equation takes 
arbitrary close to a fixed point. The meaning of the solutions associated to the fixed points can only be 
understood correctly in this way. The representation of these solutions that will be given in the following 
chapter has the only purpose of clarify the nature of the solution of such approximated equations. 
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5 Examples 


In this section we are going to apply the method described above to some specific models. We 
will start with two simple ones ( f(R ) = R n , f(R) = R + ctR n ) that are treatable also with 
the original DSA highlighting common features and differences between these two methods. 
After that, we will explore other two models (the Starobinsky and the Hu-Sawicki models) 
that are more physically interesting, but cannot be analysed with the original DSA. 


5.1 The case of l? n -gravity 

As a first check for our new formalism let us consider the first model that has been analysed 
with the original DSA [11], This model, which is also called sometimes u R n - gravity”, is 
characterized by an action in which the Ricci scalar appears as a generic power rather than 
linearly and it constitutes the simplest fourth order modification of GR. 

Putting the action in the form of Section 4.1 we have 

A = J d*xV=g [R%R n + C m ) ■ (5.1) 


Using the original DSA it was realised that for specific values of the parameter n there could be 
orbits which naturally present a transition between decelerated and accelerated expansion [11], 
This model was subsequently subjected to more detailed studies which involved cosmological 
perturbations [19, 22-24] as well as astrophysical and cosmological tests (e.g. [25-27]). The 
result of these investigations and other physical considerations points to the fact that R n - 
gravity is inconsistent with multi-scale observations and should be considered only a toy 
model. 

It is known that in terms of the original DSA, the case / oc R n is degenerate: two 
dynamical system variables (y and z) coincide. The special character of R”-gravity is present 
also in the new approach: since the term Rq appear as a factor of R n the equation for A is 
decoupled and can be excluded. Substituting in the remaining equations the expression for 
X, Y and Z 


we obtain 



4(n — 1) 

R ’ 

8 (n-2)(n — 1) 

3R2 


-— = M < 2(2 + K — M)- 

dN 1 1 ’ n- 1 


■ + I-1 

n 


= n 


3IC + (-3)M — D + 2 — 3rc 

n 


dCl 

dN 


H = 2K(K + R + 1). 


-D + l 


(5.2) 


(5.3) 


The system presents in general three invariant submanifolds: 1C = 0, D = 0 and M = 0. Table 
1 contains the standard fixed points of (5.3), together with their associated solution. The 
stability of the fixed points is shown in the case w = 0 in Table 2. 
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The solutions associated to the fixed points deserve further discussion. Since now we 
are solving the full (4.16) these solutions will be specified by a linear differential equation in 
a which is more complex than the one of the original DSA. For example for the points A, B 
and C we have 


a = Hi + a i/ 2 

a a 


H 2 sin 



+ H 3 cos 



(5.4) 


This equation can be solved exactly (but almost always implicitly) only in the case in which 
two of the constants Hi are zero. For H 2 and H 3 zero one has, for example, 


a — ao(t — to), 


(5.5) 


in the other cases the solution can only be expressed in terms of inverses of hypergeometric 
functions. Looking at the nature of equation (5.4) it is clear that for small a the power law 
behaviour is the dominant component of the solution whereas for large a the hypergeometric 
behaviour is dominant. A numerical integration of the equation (5.4) is given in Figure 1. It 


a 



Figure 1. Numerical solution of equation (5.4). The constants Hi have all been chosen to be one 
and the initial condition is a(0) = 0.01. 


is clear that the solution has sigmoid behaviour: after a power law growth of the type (5.5), 
the expansion rate starts, at first, to increase and then to decrease, to approach eventually 
a constant. This results is confirmed by the numerical check of the first derivative of the 
solution. 

The solutions associated to the points T> and £ are given by the equation 


a H 1 

— — —9-b CL 

a a z 


H 2 sin ( v3 log a) + H 3 cos ( v3 log a 


(5.6) 


As before, this equation can be solved exactly only in the case in which two of the constants 
Hi are zero. For H 2 and H 3 zero one has, for example, 


a = a 0 (t - to) 2 , 


(5.7) 
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Figure 2. Numerical solution for the energy density of (5.4). The constants H t have all been chosen 
to be one and the initial condition is a(0) = 0.01. 


a 



Figure 3. Numerical solution of equation (5.6). The constants Hi have all been chosen to be one 
and the initial condition is a(0) = 0.01. 


which is dominant only for small a. The numerical integration (Figure 3) shows that also in 
this case the solution approaches a constant a late time in a manner similar to the previous 
one. Both these solutions present two inflection points i.e. changes in the sign of the expansion 
rate of the universe. 

For T and Q instead one has, respectively, 

Cl n — 2 

- = (Hi + H 2 + H 3 ) aSP W*, (5.8) 

a 

and 

a 3(ui+i) 

- = (H 1 + H 2 + H 3 ) a -, ( 5 . 9 ) 

a 
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A* 



Figure 4. Numerical solution for the energy density of (5.6). The constants Hi have all been chosen 
to be one and the initial condition is a(0) = 0.01. 


which can be integrated exactly to give a pure power law behaviour. 

Table 1. Fixed points of f{R) = yl?” an d their associated solutions. Here ctg = Hi + H 2 + H 3 . 


Point 

Coordinates {M, K, 11} 

Scale Factor 

A 

{0,-1,0} 

(5.4) 

B 

{0,-1,-1 -3w} 

(5.4) 

C 

{n( 1 — n), 2(n — l)n — 1,0} 

(5.4) 

V 

(0,0, 2 — 3w} 

(5.6) 

£ 

{0,0,0} 

(5.6) 

T 

{jfcfe.o.o} 

(1 — 2n)(l — n) 

a = a 0 (t- t 0 ) ™- 2 

Q 

/ 3—4n+3iu n 6n 2 w-\-8n 2 — 9nw— 13n-\-3w-\-3 \ 
\ 4 n > 2 n 2 / 

2 n 

a = ao(t — t 0 ) 3(u,+1) 


It is instructive to compare the results above with the ones of [11]. It is evident that 
the fixed points we have obtained, correspond to the ones found in [11]. In particular the 
solutions associated to the fixed points are coincident when H 2 and H 3 are set to zero. This 
result shows how the original DSA would give at best an incomplete result. For example in 
(5.4), since the a^ 1 term is only relevant at small a, if A, B, C are attractors the cosmology 
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Table 2. Stability of the fixed points of f(R) = xR n in the case w = 0. Here A stays for attractor, 
R for repeller, S for saddle. 


Point 

n < |(l — a/3) 

1 (1 - y/3) < n < 0 

0 < n < 1/2 

1/2 < n < 1 

A 

S 

S 

S 

S 

B 

S 

S 

S 

s 

C 

s 

A 

S 

s 

V 

R 

R 

R 

R 

£ 

s 

S 

S 

s 

T 

A 

S 

S 

A 

Q 

S 

s 

s 

S 


Point 

1 < n < 5/4 

5/4 < n < 4/3 

4/3 < n < i (1 + V3) 

n > \ (1 + \/3) 

A 

S 

S 

S 

S 

B 

S 

S 

S 

S 

C 

s 

S 

A 

S 

V 

s 

R 

R 

R 

£ 

s 

S 

S 

S 

T 

R 

S 

S 

A 

Q 

s 

s 

s 

S 


will tend to become static after a phase of accelerated expansion. In terms of the stability, 
instead, the two methods show a complete consistency: the nature of the fixed points is 
the same as the one in [11], Specifically, comparing the ranges of stability the possibility of 
a transition between almost Friedmann (point J 7 ) and power law inflation/dark energy era 
(point Q) is present also with the new method. In fact, since both the fixed points are on 
the K = 0 invariant submanifold, we can check this result explicitly plotting this part of the 
phase space (see Figure 5). It is evident that there is, in complete agreement with [11], a set 
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of initial condition in which the transition appears. The stability analysis results in Table 2 
guarantees that this is the case also for orbits in the full phase space. 



Figure 5. Plot of a section of the invariant submanifold K = 0 for the model f(R) = X-R" • Here the 
abscissa represents the variable R and the ordinate the variable fh 


5.2 The case f(R) = R + aR n 

Let us now consider another important model for fourth order gravity: the one in which a 
generic power of the Ricci scalar is added to the Hilbert-Einstein term. Historically this is 
the most studied form of fourth order gravity because of its appearance in several quantum 
gravity calculations [28]. In cosmology the model with n = 2 has gained particular attention 
as a model of geometric inflation [4], 

The Lagrangian for this class of theories can be written in the form showed in Section 
4.1 as 

f(RoR, a) = R 0 R + aR£R n . (5.10) 

Note that for n < 0 this theory is not defined in R = 0 and in i?o = 0, but in fact we will see 
that, due to the presence of derivatives of / with respect to R in the cosmological equations 
(and therefore in the dynamical system) we will have divergences also in other intervals of n. 
Substituting the / above in the expression for X, Y and Z we obtain 


X= - + 

Y = 


(n — 


n n {AM + anG"- 1 A n R n ) ’ 


+ cm6 n-1 A ri M Tl 


Z = 


8 (n — 2 )(n — 1) 
3R2 


1 - 


+ cm6 n-1 A ri M ri 


(5.11) 
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Substituting in the general system (4.15) we obtain 


dM M [n (2 n — 3) 1C — (2n 2 + 3n + l) M + nfi + 4ro 2 — 5 n] 
dN n(n — 1) 

(.K -0 + 1 ) 

6 n_1 a(n — l)nA n_1 M n_2 ’ 


dQ 

dN 

dlC 

dN 

dA 

dN 


= 0 


3K + ( - — 3 
n 


— 0 + 2 — 3 w — 


(n — 


n( AM + an6 n 1 A n M n ) 


2K(K-R + 1), 

-2A(2 + K-M) 


(5.12) 


This system admits four invariant submainfolds (K = 0, 0 = 0, A = 0, R = 0). The last two 
invariant submanifolds can be singular, but they can still contain some fixed points. As said, 
the presence of fixed points is based on the requirement of convergence of the cosmological 
equations and the eigenvalues of the fixed points. This implies that some fixed points will 
only exist for certain values of n. 

The list of the fixed points, their associated solutions and the interval of existence of 
the fixed points are given in Table 3. Note that the points A — Q are the same of the ones 
found in the case of d? n -gravity. The presence of this type of fixed points can be understood 
thinking that nearby A = 0 the function / can be approximated with its limit for small Rq. 
For the values of n for which these fixed points exist the approximate function is RqR 71 and 
therefore the fixed points of the case / oc R^R 71 appear also in the phase space of this theory. 

On top of the points A — Q of the previous case we have some additional ones, which 
we will name Rj. These points will exist if a(n — 2) > 0 (remember that the variable A is 
defined to be non negative) and their number depends on the value of n. 

In the Ri the scale factor is described by the equation 


- = Hq + 3 H\ log a + 9H2 log 2 a, 
a 


which admits the exact solution 


a(t ) = ao exp 


^4ff 2 # 0 - H'( 
2 H 2 


tan 


\(t- t 0 )^/4H 2 H 0 - Hi 



(5.13) 


(5.14) 


For AHqH -2 — -fff > 0 this solution is monotonically growing with two inflection points at 


aresm — 


t = ip 2 = ± — 

and presents a discontinuity in 

t = t = to — 


A^HoH 2 -H‘f 


3H 2 


\/4 H 0 H 2 - H'l 


+ 2&7T, k G N. 


7r 


V 4 H 0 H 2 - H\ 


+ kir, 


ke N. 


(5.15) 


(5.16) 


For t —> t + this solution approach to zero, whereas when t —> t the solution presents a 
vertical asymptote. 
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Table 3. Fixed points of f(R) = R + aR n with their interval of existence and their associated 
solutions._ 


Point 

Coordinates {M, 1C, 11, A} 

Scale Factor 

Existence 

A 

{0,-1, 0,0} 

(5.4) 

n < 1/2 

B 

{0, —1, —1 - 3w,0} 

(5.4) 

n < 1/2 

C 

{n( 1 — n), 2 (n — l)n — 1,0,0} 

(5.4) 

n < 1 

V 

{0,0,2 - 3w, 0} 

(5.6) 

n < 1/2 

£ 

{0,0,0,0} 

(5.6) 

n < 1/2 

T 

{sfeSri.o. 0 ."} 

(1 — 2n)(l — n) 

n < 1 

a = a 0 (t-t 0 ) "- 2 

Q 

f 3— 4n+3w n 
l 4 n ’ U 

6n 2 w+8n 2 -9nw—13n-\-3w+3 

> 2iP ’ U / 

2 n 

a = ao(t — f 0 ) 3(l " +1) 

n < 1 

Ui 

|2,0,0,12 1 Va(n-2)| 

(5.13) 

a(n — 2) > 0 


In the case AH 0 H 2 — Hf < 0 the (5.13) is instead not periodic and approaches a constant. 
The expanding or contracting character depends on the sign of the quantity 


//i + \JH'i - 4 H 0 H 2 
2 H 2 


(5.17) 


P > 0 implies a growing scale factor and P < 0 a decaying one. Plots of the solution (5.13) 
can be found in Figure 6. 

In the case AHqH 2 — Hf > 0 the solution (5.14) presents features which are physically 
very interesting. For times close to t + the solution grows exponentially, then around t\ it 
changes into a decelerated expansion with non-constant deceleration factor. After a coasting 
phase around t^, the decelerated expansion is followed by a new accelerated expansion phase. 
In a finite time (at t = t~), however, the solution becomes singular in the sense that a and 
its derivatives as well as the Ricci scalar diverges at this specific time. This is a well known 
property of /(i?)-gravity 7 [31, 32], but in this context one is able to appreciate both this 
drawback of the theory and its potential as a model that unify inflation and dark energy. 

Note that for H\ 2 = 0 the solution (5.14) reduce to the standard de Sitter solution. 
This fact on one hand connects the points 77; to the de Sitter fixed point of the original DSA. 
On the other hand gives a hint of the true meaning of de Sitter solutions in the framework of 
f(R)- gravity. 

' These models in fact present also other types of singularities, like the “weak singularities” in [33] or the 
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0.5 1.0 1.5 2.0 2.5 3.0 


3.5 


(a) Plot of a period of (5.13) in the case AHqHo — H\ > 0. Here 
Hi = 1, H 2 = 3, H 3 = 2 and ao = 1. 


(b) Plot of (5.13) in the case AH 0 H 2 — Hi < 0. The 
constants Hi have all been chosen so that P > 0 and 
ao = 1 


(c) Plot of (5.13) in the case 4H 0 H 2 - Hi < 0. The 
constants Hi have all been chosen so that P < 0 and 
a 0 = 1 

Figure 6. Plots of (5.13) illustrating the different behaviour that this solution can represent. 




ones found in [34]. 
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Figure 7. Numerical solution for the energy density of (5.13) in the case 4:HqH 2 — Hf > 0. Here 
Hi = 1, H 2 = 3, H 3 = 2 and cio = 1. 


The stability of the hxed points can be calculated as in the previous case using the 
Hartman Grobman theorem, and it is illustrated in Table 4. It is tempting to attempt a 
general derivation for the stability of points Hi, however numerical inspection shows that 
such analysis might be unreliable. The calculations show that when points Hi exist there are 
only two options for their stability: they can be either a saddle or an attractor (or their focus 
counterpart). In Table 4 we report some samples of the stability of points Hi for different 
values of the parameters a and n. Examples of the phase space for these theories in the form 
of the invariant submanifold R, A is given in Figures 8. Since this model can be treated 
also with the original DSA it is useful to make a comparison between the results we obtained 
above and the ones in [14], Differently from the case of R n gravity the phase space obtained 
by the two methods is not the same. In particular, the original DSA returns a phase space 
with many more fixed points. The origin of this difference is probably to be attributed to the 
choice of variables of the original DSA. For the common points one can compare the results on 
the stability and it is easy to verify complete consistency. For example, the de Sitter solution 
that in [14] is associated with £* has a stability that coincide with the point H of the present 
analysis. 

In the case n = 2 the new DSA returns a phase space with no finite hxed points. It is 
known that the case n = 2 in the theory (5.10) presents significant physical differences with 
respect to the other model of this class [28], and it is to be expected that the phase space will 
reflect these differences. The fact that the phase space does not present a point of type H 
does not necessarily imply that the model does not have de Sitter solutions (we know in fact 
that they are present [30]). A careful analysis of the equations shows that the hxed point in 
this case is asymptotic (A —> oo) and it is therefore excluded by the present analysis. 

It is interesting to note that one of the conclusions in [31] is that to avoid the singularity 
one either has to recur to special initial conditions or to add additional curvature invariants. 
This result seems consistent with our findings. In order to avoid the singularity one either has 
to control the initial conditions (by setting to zero some of the constant Hi), the values of the 
parameters like in Figure 5.2, or hope that adding additional curvature invariants the hxed 
points H will become irrelevant (in the same way of what happens with the case n = 2 above). 
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Table 4. Stability of the fixed points of f(R) = R + aR n in the case w = 0. Here n\ is the smallest 
real solution of the equation 256n 3 — 608n 2 + 417?r — 81 = 0 , A stays for attractor, R for repeller, S 
for saddle, Fg for saddle focus. The fixed points appear in the table only if they exist in at least one 
of the intervals of the parameter n indicated. 


Point 

n < \ (l — \/3) 

l (1 - V3) < n < 0 

0 < n < n\ 

n\ < n < 1/2 

A 

S 

S 

S 

S 

B 

F S 

F S 

f 5 

Fs 

C 

S 

S 

s 

S 

V 

S 

S 

s 

S 

£ 

s 

s 

s 

s 

T 

A 

s 

s 

s 

Q 

F S 

F S 

F 5 

s 


Point 

1/2 < n < 1 

n > 1 

C 

S 

NA 

T 

S 

NA 

Q 

F S 

NA 


Table 5. Stability of the fixed points 77; of f(R) = R + aR n in the case w = 0. Here A stays for 
attractor, S for saddle, for saddle focus and NA represents the absence of fixed points. In the 
cases considered there is only one fixed point 77. 



n = —2 

n = —1 

n = 3/2 

n = 5/2 

a = — 2 

S 

S 

NA 

NA 

a = — 1 

S 

S 

NA 

NA 

a = 1 

A 

A 

NA 

NA 

a = 2 

NA 

NA 

S 

S 
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(a) Plot of the R > 0 section of the phase space in 
the case n = l/4a = —10. 



(b) The case n = 3 a = 1/16. 


Figure 8. Samples of the invariant submanifold Q = 0, K = 0 for the theory f(R) = R + aR n . 
The values of the parameters have been chosen to give the best graphical representation of the phase 
space. 


The issue is that the feature of the inflation dark energy connection seems to be inextricably 
tied to the approach to the singularity. Thus, in order to generate cosmic acceleration without 
incurring in the singularity requires that other dynamical mechanisms/forms of the function 
/ have to be found. 

5.3 The Starobinsky model. 

The Starobinsky model is one of the most important class of models of Dark Energy based on 
fourth order gravitation [35]. The basic idea is to construct, via a function of the Ricci scalar, 
an effective cosmological constant which is relevant in curved spacetimes, but approaches zero 
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in the case of flat spacetime. The function / of this model is given by 


/ (Ro, R) — R + Ai?o 



(5.18) 


with n and A positive and Rq of the order of the inverse of the present value of the cosmological 
constant. The parameter A and the value of n are related by an algebraic equation. In [35] 
considerations on the Solar System constraints and cosmological linear perturbation theory 
were used to find that one should have n >2 and A > 0.94. 

Following Section 4.1 we can write the Lagrangian (5.18) for this class of models as 


f(R) = RqR + a 


(i + p 2 r%r 2 ) n -i 


(5.19) 


where the parameters a,/3 are related to the ones in (5.18) by the relations A = a \(3\ Ro and 

Ro = mR 0 r 1 . 

The functions X, Y, Z are 


X = 

Y = 

Z = 


a + 12a/3AM[3AM(l + 2 n) — an] 6AM — a 


6A 


(1 + 6 2 A 2 /3M 2 ) n+1 - 12aA/?nM 


3 2 4 ct/3nA (36A 2 /3(2n + 1)M 2 - l) 


6A 


(1 + 6 2 /3A 2 M 2 ) (1 + 6 2 A 2 /3M 2 ) n+1 - 12aA/3nM 
2 8 3 3 a/3 2 n(n + 1)A 5 M (l2A 2 /3(2n + 1)M 2 - l) 


(1 + 6 2 /?A 2 M 2 ) 2 (1 + 6 2 A 2 /3M 2 ) n+1 - 12aA/3nM 


(5.20) 


and the dynamical system (4.15) becomes 


dM. 

dN 


dQ 

dN 


dK 

dN 

dA 

dN 


26 2 nA 2 (2n6 2 A 2 M 2 + 6 2 A 2 M 2 - 1) ^ 1 + 2 M + + 3) + M] 

6 4 /3 2 M 3 A 4 [2n(3M + fi - K(4n + 3) + 4n(M - 2) - 5) + M]} 

[6A(K - Q + 1) - a] (6 2 A 2 /3M 2 + l)” +2 
26 2 a/3nA 2 [36A 2 ^(2n + 1)M 2 - 1] ’ 

afI{12/3AM[3AM(2n + 1) - an] + 1} 


6A 


(1 + 6 2 A 2 /3M 2 ) 


n+1 


12aA/3nM 


(5.21) 


n{A[6(2 + 3w ) - 18K + 12 M + 611] + a} 
6A 


2K(K — M + 1), 


-2A(2 + K-R). 


The system presents three invariant submanifolds: K = 0, 12 = 0, A = 0, although the last 
submanifold is singular. Note that, differently form the previous case, none of the fixed points 
A — G is present here. This can be explained looking at the limit of the action for Ro —> 0. In 
this case in fact the action reduces to RqR rather than RqRJ 1 and therefore the fixed points 
of i?"-gravity are not present. 
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Even if this system is valid for any value of the parameters, it is not possible to find its 
fixed points analytically 8 : this task would entail the resolution of an algebraic equation of 
order n for which no general analytical solution is known. We therefore refer to the bounds 
in [35] and we set from now on n = 2 and f3, a = 1. In this case the system admits three 
fixed points with real coordinates (see Table 6). One of them is on the singular submanifold 
and presents singular eigenvalues so that it does not appear in Table 6. The other two (77i 
and 7 ^ 2 ) are both associated to a solution of the type (5.13). 

The stability analysis reveals that the character of the points 77* is in general different. 
In our particular case, only one of these point is an attractor whereas the other is unstable. 
Since the phase space contains invariant submanifolds, we can conclude that for these values of 
the parameters there is only a specific set of initial conditions which lead to the attractor 772- 
Therefore, also in the case of the Starobinsky model there is the possibility that the cosmology 
will evolve towards a singularity, but this occurrence depends strictly on the choice of the 
initial conditions. In orbits which do not approach 77, the fate of the cosmological models 
depends on the presence of asymptotic attractors. 

Table 6. Fixed points for the Starobinsky model and their stability in the case n = 2 and a = j3 = 1. 
Here A 0 is the only positive real solution of the equation (6A — 1) (144A 2 + l) 4 + 144A 2 (432A 2 + 4) + 
1 = 0 different form 1/12, S represents a saddle and is an attractive focus. 


Point 

Coordinates {M,IK, 17, A} 

Scale Factor 

Stability 

77 1 

Ho 

{2,0,0,^} 

{2,0,0, A 0 } 

(5.13) 

(5.13) 

S 

Fa 


A plot of the section of the invariant submanifold IK = 0,12 = 0 which contains 7~L\ and 
7^2 is given in Figure 9. 

5.4 The Hu-Sawicki model 

As a last example we consider a model for geometric DE proposed by Hu and Sawicki [36]. 
The model was designed to be able to reproduce cosmic acceleration without the explicit intro¬ 
duction of a cosmological constant and, at the same time, to be compatible with cosmological 
and Solar System tests. 

The action for the Hu-Sawicki model can be written as (4.1) setting 

rv Tt n 

f(R 0 ,R) = R 0 R - r? -, (5.22) 

’ ' 1 + /37^7K v 

where the parameter n is chosen to be positive, a > 0 and f3 > 0. With this choice of / the 

8 Of course the calculation could be done numerically, but we do not perform such task here. 
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0.20 



Figure 9. The section of the invariant submanifold K = 0, Q = 0 of the phase space of the Starobinsky 
model which contains Hi and H 2 in the case n = 2, Rq > 0 and (3,a = 1. 


functions X, Y, Z are given by 

v _ 6/3AM - a | a 6 2 A 2 M 2 [/36 n A n M n (l + n) - cm6 n - 1 A n - 1 M n - 1 + l] 
_ 6A/3 f 6AM + 6 n A n M n (12A/3M - an) + /3 2 6 2n+1 A 2n+1 M 2n+1 ’ 

_ 4n [12AM — 6 n A n M n (a — 12A/3M + cm)] 

_ 6AM 2 [6AM + 6 n A n M n (12A/3M - an) + /3 2 6 2n+1 A 2n+1 M 2n+1 ] 

8 n 

~~ M (/36 n A n R n + 1) ’ 

2 3 n6 2 AM [/36 n A”M n (n + 1) - a6 n - 2 A n " 1 M n - 1 (n 2 + 2 + 3n) + l] 

A 3M 2 [6AM + 6 n A n M n (12A/3M - cm) + /3 2 6 2n+1 A 2n+1 M 2n+1 ] 

16n 2 16n (n + 1) 

+ M^G^A^Tl? ~~ M2^36^M^Tl)' 
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Substituting in (4.15), we obtain the dynamical system 


c® 

dN 


d£l 

dN 

dlC 

dN 

dA 

dN 


- - - -k {cm(3K — Q + 5) + M[a — 12A/3(K — fi + 1)1 + 2cm 3 (K — M + 2) 

anyn + 1 + 

+cm 2 (5K — 5M — 14 + 9) — 24nA/3M(I€ — 11 + 1)} 


+ 


2nM{[24A/3nM - aK(n 2 - 1)](K - ft + 1) + a(n - 1) 2 M} 


a/36 n (n — l)(n + l) 3 A n M n — a (n 2 — 1)^ 


+ 


/36 n A n M n+2 [a - 6A/3(1 + Q + K)] 6 1-n A 

an(n + 1) 

aMfl [(/3(n + 1) — na)6 n nA n W 1 + 1] 


l-n a l-nin>2-n 


6/3AM (/3 6 n A n M n + 1) - a/3n6 n A n R n 


an(n — 1) 

-n(n-3K + 


(K-n + i), 


+ 3w — 2 + 


a 


6A/3 


2K(K — M + 1), 


-2A(2 + K-R), 


(5.24) 


This system can present four invariant submanifolds (K = 0, H = 0, A = 0 and M = 0), 
but the last one exists only for 0 < n < 2. The A = 0 submanifold can be singular. This 
happens for n > 1. Differently from the Starobinsky model, the phase space in this case 
contains different fixed points depending on the value of the parameter n. In particular, the 
conditions /3 > 0 and A > 0 limits strongly the number of fixed point belonging to the class 
%. The list of fixed points for n / 1 is given in Table 7. The case n = 1 requires a special 
treatment and will be covered in a separate subsection. 

As usual, the conditions of existence for the fixed points in Table 7 have been determined 
asking that the dynamical system should be at least of class C{ 2) in the fixed point. In 
one case (the point C ) there are values of n for which the Jacobian is divergent, but the 
eigenvalues converge. We refer as interval of existence of the point C the one of convergence 
of the eigenvalues. Within the interval of existence of the fixed points the stability can be 
determined with the standard methods. The results are given in Table 8. As in the case of 
f(R) = R + aR n the stability of 77,; needs to be calculated case by case, but it can only be 
a saddle or an attractor of their focus counterparts. An example of the stability of points Hi 
for different combinations of the values of a, /3 and n in the case of dust (w = 0) is given in 
Table 8 

In this case the points C and 77 are the only possible finite attractors for the cosmology. 
The first point, however, only exists for 0 < n < 1. Therefore for n > 1, we are in the same 
situation of the Starobinsky model: only a careful choice of the initial conditions could avoid 
the singularity of solution (5.13). 


5.4.1 The case n = 1. 

It is interesting to note that the dynamical system (5.24) is not defined for n = 1, which 
means that the dynamical system in this case needs to be obtained re-deriving X, Y, Z and 
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Table 7. Fixed points of the Hu-Sawicki model for n ^ 1 and their associated solutions. Here Ao 
represents the positive solutions of the equation given in the last line of the Table. 


Point 

Coordinates {M, 1C, 17, A} 

Scale Factor 

Existence 

A 

{0,-1, 0,0} 

(5.4) 

0 < n < 1 

B 

{0,-1,-1-310,0} 

(5.4) 

0 < n < 1 

C 

{n( 1 — n), 2(n — l)n — 1,0,0} 

(5.4) 

0 < n < 1 

V 

{0,0,2 - 3w, 0} 

(5.6) 

0 < n < 1 

£ 

{0,0,0,0} 

(5.6) 

0 < n < 1/2 

T 

O 

o' 

cT 

1 1 

(1 —2n)(l- 

-n) 

a = a 0 (t-t 0 ) "~ 2 

0 < n < 1 

'Hi 

{2,0,0, A 0 } 

(5.13) 


12- ,l A" 

" n {l2A+2 6n+1 /3 2 27”A 3n (6A/3—a)+/3144 n A 2n [36A/3+a(n—4)]+12”A” 

[36A/3+a(n—2)]} _ n 


Oin[(312 n (n-\-l)A n — n+1] 



the equations (4.15). This procedure yields 


2(3/1 AM + 1)’ 


M(6/3AM + 1) M(3/3AM + 1) ’ 

32 | 16 16 

~M 2 (6/3AM + 1) + M 2 (6/3AM + l) 2 + M 2 (3/3AM + 1)’ 


(5.25) 
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Table 8. Stability of the fixed points of the Hu-Sawicki model in the case n / 1 and a / 1. Here A 
stays for attractor, S for saddle and for attractive focus. The value of the A coordinate of R is 
approximated. 


Point 0 < n < 1/2 

1/2 < n < 1 

n > 1 

A S 

S 

NA 






(n, a, (3) 

Coordinates of Hi 

Stability 

B 

S 

S 

NA 




C 

Fa 

Fa 

NA 

(2, 3,1/2) 

{2,0,0,0.07} 

{2,0,0,0.97} 

S 

Fa 

V 

S 

S 

NA 

(3,4,7/10) 

{2,0,0,0.89} 

{2,0,0,0.95} 

Fa 

S 

£ 

S 

NA 

NA 




T 

S 

S 

NA 

(4, 3, 3/5) 

{2,0,0,0.10} 

{2,0,0,0.83} 

S 

Fa 

Q 

S 

S 

NA 





and, therefore, 


{ aK(1 + 30/3AM) " (K + 1)(6/3AR +1)3 
+a[18/3AM(2/3AM 2 - R + 3) + 1]} 
n(6/3AM + 1) [1 - a + 12/3 AM(1 + 3/3AM)] 

+ 12a/3A ’ 

W = „-(8<s£» + 1 )» {( ® - 2R - ^ + 2 )I“ - 1 - 12A « 3A '® + D] ( 5 26 > 

+6aA/3M 2 } — n 2 , 

^ = 2K(K-R + 1), 

fik 

-— = -2A(2 + K-R). 


Differently from the system (5.24), (5.26) does not present the invariant submanifold 
M = 0, and the submanifold A = 0 is singular. The properties of the phase space depends on 
the values of a. 

In particular, for a/1 the phase space contains three fixed points of the type 7~L and an 
additional one T~Lq = {M = 2, K = 0, D = —(4 + 3w), A = — prg}. However, all of these points 
have A < 0 for a and /3 positive and have to be discarded. Therefore, this version of the Hu- 
Sawicki model presents a finite phase space structure analogous to the one of f(R) = R+aR 2 . 
The difference is that the Hu-Sawicki model dose not necessarily have asymptotic fixed points 
of the type R. 
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For a = 1, the situation is very different. In this case the function / reduces to 


f(R,Ro) 


/3R 0 R 
1 + /3RqR 


(5.27) 


i.e. the Hilbert Einstein term is eliminated. The phase space for this case contains the 
points A (which for n = 1 coincides with C ) to £. Three fixed points of the type R exist, 
but two of them have A < 0 for /3 positive and have to be discarded. The third one has 
coordinates {M = 2,1C = 0,11 = 0, A = 0}. In addition, two new fixed points with coordinates 
X = {M = 4,K = 3,H = 0,A = 0} and C = {R = £(5 - 3iu),K = 0,0 = -§(1 + w),A = 0} 
are present. 

The modified Raychaudhuri equation (4.16) returns the same solutions for the fixed 
points A-£. The solution associated to X is (5.4) whereas the one associated to £ is 

a = ao(t — to) 3(1+ “ ) ■ (5.28) 


The stability of the fixed points in this case can be calculated in the standard way and 

with the exception of R all appear to be unstable (see Table 9). Point R has instead one zero 

eigenvalue and the analysis of its stability requires the use of the Center Manifold Theorem 
(CMT) [37]. To apply this theorem the first step is to write the system (5.26) in the form 

CA + F(A,Y) (5.29) 

PY + G(A,Y) (5.30) 


dA 

dN 

dY 

dN 


where Y = {1C, M, H}, M = 3(2 — M) — 12/3A, C corresponds to the linear pat of the equation 
for A, the vector P to the linear part of the equation of Y and F and the vector G represent 
the non-linear part of the equation for A and Y. The CMT tells us that the behaviour of the 
fixed points is determined by the solution h of the equation 


dA 

dN 


CA + F(A, h(A)) 


(5.31) 


where the vector function h(A) is given by 


h'(A) [CA + F( A, h( A)] - [P/i(A) + G(A, h{ A))] = 0 (5.32) 


Approximating the function h(A) with its Taylor series one finds that 

h(A) = aA 2 + ... 


(5.33) 


where a is a vector that depends on the parameter (3. Since the first non-zero term in the 
previous expression is quadratic, (5.31) implies that point R has the stability of a saddle-node. 
In Figure 10 we plot an example of the invariant submanifold 1C = 0, H = 0 corresponding to 
this case. 


6 Conclusions 

In this paper we have presented a new approach to analyse the finite phase space of the 
cosmology of f(R)- gravity. The new method can be applied to any form of the function / 


- 28 - 




Table 9. Stability of the fixed points of f the Hu-Sawicki model in the case n = 1 and a = 1. Here 
A stays for attractor, S for saddle, for attractive focus and SN for non hyperbolic saddle-node. 

Point Coordinates {M, K, fl, A} Scale factor Stability 


A 

{0,-1,0,0} 

(5.4) 

S 

B 

{0,-1,-1 -3w,0} 

(5.4) 

S 

V 

{0,0,2 - 3w,0} 

(5.6) 

S 

£ 

{0,0, 0,0} 

(5.6) 

S 

n 

{2, 0,0,0} 

(5.6) 

SN 

X 

{4, 3, 0,0} 

(5.13) 

S 

£ 

{i(5-3u;),0,-§(l + u7),0} 

a = ao{t — to) 3(1+u,) 

s 



Figure 10. The section of the invariant submanifold K = 0, fl = 0 in the case n = 1 , a = 1 , /3 = T 
for the Hu-Sawicki model. Note that on the A > 0 part of this invariant submanifold TL appears an 
attractor, but has a saddle character in the unphysical A < 0 of the phase space (not represented 
here). 


which is analytical in R without the need of cumbersome inversions. The phase space obtained 
is more regular than the one obtained with the original dynamical systems approach, although 
it is not possible to eliminate singularities in a complete way as they are due to the essential 
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form of the cosmological equations. The new method also naturally excludes the fixed points 
which represent states incompatible with the definition of the dynamical system variables. 
Therefore the new DSA returns a phase space which matches in a closer way the actual 
evolution of the cosmological equations. 

Using the idea of higher order cosmological parameters (like “jolt” and “snap”), the 
new DSA is able to associate to the fixed points full solutions of the cosmological equations 
(in the sense of solutions with four integration constants). Among the fixed points found 
in our analysis, the points labeled Hi are surely the most interesting. They correspond 
to the dominance of the curvature terms and their number depends on the value of key 
parameters appearing in the function /. The scale factor in H is a transcendental function 
with the remarkable property to be able to combine an initial exponential expansion a phase 
of decelerated expansion and a final accelerated expansion phase. Its existence on one hand 
describes clearly the role of /(i?)-gravity as model for double inflation or a unified model for 
inflation, standard cosmology and dark energy. On the other, however, it confirms clearly 
that f(R) cosmologies with these properties can run into finite time singularities. This is 
indicated by the fact that in many of the cases we have analysed, the solution H is the only 
attractive fixed point of the phase space (although not a global attractor), ft is easy to 
prove with our method, and in accordance to the results in literature (see e.g. [31, 33]) that 
adding a special set of additional invariants one might avoid this scenario, but in general such 
avoidance requires fine-tuning. 

ft is important to be careful in considering the nature of the points H. Their associated 
solution (5.13) contains integration constant which can take any value and it is in general 
related to the constants appearing in /. This means that the fact that a theory posses an 
attractor H does not mean that the full behaviour (5.13) is necessarily realised: one could 
have, for example, that for specific values of the coupling constants and of the Hi vanishes. 
If H 2 and H% are zero then H represents a standard de Sitter solution. This depends on the 
value of the coupling constant of the model as well as the initial condition for the model. 
However it is evident that the nature of (5.13) has repercussions on the understanding of the 
actual meaning and role of the de Sitter solutions in the context of f(R)~ gravity. 

The DSA proposed has first been verified on the simple case of f{R) = R n and we found 
a complete agreement in terms of fixed points and stability with the original DSA in [11], 
However, the solutions associated to the fixed points present significative differences. In fact, 
even if the solutions obtained with the new method reduce to the ones of the old DSA setting 
to zero two of the four integration constants, it appears clear that their (time-)asymptotic 
behaviour is different. This means that, for example, when the fixed points are attractors, the 
parts of the solution excluded in the original DSA can be dominant. In this respect, therefore, 
the original DSA returns incomplete information on the cosmology of f(R)~ gravity. 

As another test of the new DSA, we have considered another form of / that was analysed 
with the original dynamical system method: f(R) = R + aR n . In this case, it turns out that 
the only difference between the two treatment is in the number of fixed points: the new method 
excludes fixed points that give conditions inconsistent with the definition of the variables. The 
remaining fixed points, however, present a stability that is completely equivalent. The theory 
f(R) = R + aR n is also the simplest theory that shows points of the type R in the finite 
phase space. An interesting exception in this respect is the case n = 2 and it is natural to 
expect that this difference is related somehow to the special properties of this model. 

As last steps we have applied the new DSA to two models which were not analysable 
with the original method, but at the same time constitute important theoretical models for 
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inflation and/or dark energy: the Starobinsky and the Hu-Sawicki models. 

The Starobinsky model, unfortunately, cannot be treated in general due to the complex¬ 
ity of the algebraic equations needed to determine the fixed points. We have therefore limited 
our analysis to the case in which the value of the parameters are chosen to be compatible 
with Solar System and cosmological perturbations constraints given in [35]. In this specific 
case, we obtain a phase space in which only two fixed points of the type PL appear. Since 
only one of the points is an attractor we can conclude that the scenario of solution (5.14) is 
possible in this model, but it is not achievable for general initial conditions. 

The case of the Hu-Sawicki model is more involved. Only for specific intervals of the 
parameter n the phase space admits fixed points different form the type PL. The situation 
is complicated by the fact that the points belong to singular submanifold. The case n = 1 
has to be treated separately and reserves a number of surprises. For example, in the case 
n = 1, a / 1 the phase space has no finite fixed points, much in the same way of the case 
f(R) = R + aR 2 . The two models, however, differ in the asymptotic structure of the phase 
space. In the case n = 1, a = 1 we found, instead, that only the (unique) point PL is always 
effectively an attractor. 

The new DSA seems therefore to be very efficient in uncovering the features of interesting 
f(R) cosmologies. However, as all methodologies, the new DSA presents also a series of 
drawbacks. For example, one would like to be able to consider in an easier way the GR-like 
states for the cosmology to be able to find (if they exist) "mimicking behaviours" of these 
theories, not dissimilar to the isotropization mechanism already found in scalar tensor gravity 
[38, 39]. In addition, the impossibility to define a set of compact variables and the fact that 
our results point clearly to the presence of asymptotic fixed points makes the present analysis 
incomplete. The resolution of these issues in the context of the new DSA will be the focus of 
future studies. 
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